<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dgels.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="DGELS.1"></a><a href="dgels.f.html#DGELS.1">DGELS</a>( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
     $                  INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK driver routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          TRANS
      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DGELS.19"></a><a href="dgels.f.html#DGELS.1">DGELS</a> solves overdetermined or underdetermined real linear systems
</span><span class="comment">*</span><span class="comment">  involving an M-by-N matrix A, or its transpose, using a QR or LQ
</span><span class="comment">*</span><span class="comment">  factorization of A.  It is assumed that A has full rank.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The following options are provided:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  1. If TRANS = 'N' and m &gt;= n:  find the least squares solution of
</span><span class="comment">*</span><span class="comment">     an overdetermined system, i.e., solve the least squares problem
</span><span class="comment">*</span><span class="comment">                  minimize || B - A*X ||.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  2. If TRANS = 'N' and m &lt; n:  find the minimum norm solution of
</span><span class="comment">*</span><span class="comment">     an underdetermined system A * X = B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  3. If TRANS = 'T' and m &gt;= n:  find the minimum norm solution of
</span><span class="comment">*</span><span class="comment">     an undetermined system A**T * X = B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  4. If TRANS = 'T' and m &lt; n:  find the least squares solution of
</span><span class="comment">*</span><span class="comment">     an overdetermined system, i.e., solve the least squares problem
</span><span class="comment">*</span><span class="comment">                  minimize || B - A**T * X ||.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Several right hand side vectors b and solution vectors x can be
</span><span class="comment">*</span><span class="comment">  handled in a single call; they are stored as the columns of the
</span><span class="comment">*</span><span class="comment">  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
</span><span class="comment">*</span><span class="comment">  matrix X.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TRANS   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'N': the linear system involves A;
</span><span class="comment">*</span><span class="comment">          = 'T': the linear system involves A**T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of rows of the matrix A.  M &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of the matrix A.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NRHS    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of right hand sides, i.e., the number of
</span><span class="comment">*</span><span class="comment">          columns of the matrices B and X. NRHS &gt;=0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
</span><span class="comment">*</span><span class="comment">          On entry, the M-by-N matrix A.
</span><span class="comment">*</span><span class="comment">          On exit,
</span><span class="comment">*</span><span class="comment">            if M &gt;= N, A is overwritten by details of its QR
</span><span class="comment">*</span><span class="comment">                       factorization as returned by <a name="DGEQRF.65"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>;
</span><span class="comment">*</span><span class="comment">            if M &lt;  N, A is overwritten by details of its LQ
</span><span class="comment">*</span><span class="comment">                       factorization as returned by <a name="DGELQF.67"></a><a href="dgelqf.f.html#DGELQF.1">DGELQF</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A.  LDA &gt;= max(1,M).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
</span><span class="comment">*</span><span class="comment">          On entry, the matrix B of right hand side vectors, stored
</span><span class="comment">*</span><span class="comment">          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
</span><span class="comment">*</span><span class="comment">          if TRANS = 'T'.
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, B is overwritten by the solution
</span><span class="comment">*</span><span class="comment">          vectors, stored columnwise:
</span><span class="comment">*</span><span class="comment">          if TRANS = 'N' and m &gt;= n, rows 1 to n of B contain the least
</span><span class="comment">*</span><span class="comment">          squares solution vectors; the residual sum of squares for the
</span><span class="comment">*</span><span class="comment">          solution in each column is given by the sum of squares of
</span><span class="comment">*</span><span class="comment">          elements N+1 to M in that column;
</span><span class="comment">*</span><span class="comment">          if TRANS = 'N' and m &lt; n, rows 1 to N of B contain the
</span><span class="comment">*</span><span class="comment">          minimum norm solution vectors;
</span><span class="comment">*</span><span class="comment">          if TRANS = 'T' and m &gt;= n, rows 1 to M of B contain the
</span><span class="comment">*</span><span class="comment">          minimum norm solution vectors;
</span><span class="comment">*</span><span class="comment">          if TRANS = 'T' and m &lt; n, rows 1 to M of B contain the
</span><span class="comment">*</span><span class="comment">          least squares solution vectors; the residual sum of squares
</span><span class="comment">*</span><span class="comment">          for the solution in each column is given by the sum of
</span><span class="comment">*</span><span class="comment">          squares of elements M+1 to N in that column.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDB     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array B. LDB &gt;= MAX(1,M,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK.
</span><span class="comment">*</span><span class="comment">          LWORK &gt;= max( 1, MN + max( MN, NRHS ) ).
</span><span class="comment">*</span><span class="comment">          For optimal performance,
</span><span class="comment">*</span><span class="comment">          LWORK &gt;= max( 1, MN + max( MN, NRHS )*NB ).
</span><span class="comment">*</span><span class="comment">          where MN = min(M,N) and NB is the optimum block size.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.107"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  if INFO =  i, the i-th diagonal element of the
</span><span class="comment">*</span><span class="comment">                triangular factor of A is zero, so that A does not have
</span><span class="comment">*</span><span class="comment">                full rank; the least squares solution could not be
</span><span class="comment">*</span><span class="comment">                computed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, TPSD
      INTEGER            BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMLNUM
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      DOUBLE PRECISION   RWORK( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.132"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.133"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      DOUBLE PRECISION   <a name="DLAMCH.134"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.134"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
      EXTERNAL           <a name="LSAME.135"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.135"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>, <a name="DLABAD.135"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, <a name="DLAMCH.135"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.135"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DGELQF.138"></a><a href="dgelqf.f.html#DGELQF.1">DGELQF</a>, <a name="DGEQRF.138"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>, <a name="DLASCL.138"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASET.138"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, <a name="DORMLQ.138"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</a>, <a name="DORMQR.138"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>,
     $                   <a name="DTRTRS.139"></a><a href="dtrtrs.f.html#DTRTRS.1">DTRTRS</a>, <a name="XERBLA.139"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          DBLE, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input arguments.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      MN = MIN( M, N )
      LQUERY = ( LWORK.EQ.-1 )
      IF( .NOT.( <a name="LSAME.151"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> ) .OR. <a name="LSAME.151"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'T'</span> ) ) ) THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -6
      ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
         INFO = -8
      ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
     $          THEN
         INFO = -10
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Figure out optimal block size
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
<span class="comment">*</span><span class="comment">
</span>         TPSD = .TRUE.
         IF( <a name="LSAME.173"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> ) )
     $      TPSD = .FALSE.
<span class="comment">*</span><span class="comment">
</span>         IF( M.GE.N ) THEN
            NB = <a name="ILAENV.177"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DGEQRF.177"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            IF( TPSD ) THEN
               NB = MAX( NB, <a name="ILAENV.179"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DORMQR.179"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>'</span>, <span class="string">'LN'</span>, M, NRHS, N,
     $              -1 ) )
            ELSE
               NB = MAX( NB, <a name="ILAENV.182"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DORMQR.182"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>'</span>, <span class="string">'LT'</span>, M, NRHS, N,
     $              -1 ) )
            END IF
         ELSE
            NB = <a name="ILAENV.186"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DGELQF.186"></a><a href="dgelqf.f.html#DGELQF.1">DGELQF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            IF( TPSD ) THEN
               NB = MAX( NB, <a name="ILAENV.188"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DORMLQ.188"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</a>'</span>, <span class="string">'LT'</span>, N, NRHS, M,
     $              -1 ) )
            ELSE
               NB = MAX( NB, <a name="ILAENV.191"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="DORMLQ.191"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</a>'</span>, <span class="string">'LN'</span>, N, NRHS, M,
     $              -1 ) )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
         WORK( 1 ) = DBLE( WSIZE )
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.202"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DGELS.202"></a><a href="dgels.f.html#DGELS.1">DGELS</a> '</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( MIN( M, N, NRHS ).EQ.0 ) THEN
         CALL <a name="DLASET.211"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine parameters
</span><span class="comment">*</span><span class="comment">
</span>      SMLNUM = <a name="DLAMCH.217"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) / <a name="DLAMCH.217"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      BIGNUM = ONE / SMLNUM
      CALL <a name="DLABAD.219"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( SMLNUM, BIGNUM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale A, B if max element outside range [SMLNUM,BIGNUM]
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="DLANGE.223"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>( <span class="string">'M'</span>, M, N, A, LDA, RWORK )
      IASCL = 0
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm up to SMLNUM
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASCL.229"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
         IASCL = 1
      ELSE IF( ANRM.GT.BIGNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm down to BIGNUM
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASCL.235"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
         IASCL = 2
      ELSE IF( ANRM.EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Matrix all zero. Return zero solution.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASET.241"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'F'</span>, MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
         GO TO 50
      END IF
<span class="comment">*</span><span class="comment">
</span>      BROW = M
      IF( TPSD )
     $   BROW = N
      BNRM = <a name="DLANGE.248"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>( <span class="string">'M'</span>, BROW, NRHS, B, LDB, RWORK )
      IBSCL = 0
      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm up to SMLNUM
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASCL.254"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
     $                INFO )
         IBSCL = 1
      ELSE IF( BNRM.GT.BIGNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm down to BIGNUM
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASCL.261"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
     $                INFO )
         IBSCL = 2
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( M.GE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        compute QR factorization of A
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DGEQRF.270"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
     $                INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        workspace at least N, optimally N*NB
</span><span class="comment">*</span><span class="comment">
</span>         IF( .NOT.TPSD ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Least-Squares Problem min || A * X - B ||
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DORMQR.281"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>( <span class="string">'Left'</span>, <span class="string">'Transpose'</span>, M, NRHS, N, A, LDA,
     $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           workspace at least NRHS, optimally NRHS*NB
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DTRTRS.289"></a><a href="dtrtrs.f.html#DTRTRS.1">DTRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non-unit'</span>, N, NRHS,
     $                   A, LDA, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>            IF( INFO.GT.0 ) THEN
               RETURN
            END IF
<span class="comment">*</span><span class="comment">
</span>            SCLLEN = N
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Overdetermined system of equations A' * X = B
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DTRTRS.304"></a><a href="dtrtrs.f.html#DTRTRS.1">DTRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'Transpose'</span>, <span class="string">'Non-unit'</span>, N, NRHS,
     $                   A, LDA, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>            IF( INFO.GT.0 ) THEN
               RETURN
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(N+1:M,1:NRHS) = ZERO
</span><span class="comment">*</span><span class="comment">
</span>            DO 20 J = 1, NRHS
               DO 10 I = N + 1, M
                  B( I, J ) = ZERO
   10          CONTINUE
   20       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DORMQR.321"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>( <span class="string">'Left'</span>, <span class="string">'No transpose'</span>, M, NRHS, N, A, LDA,
     $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           workspace at least NRHS, optimally NRHS*NB
</span><span class="comment">*</span><span class="comment">
</span>            SCLLEN = M
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute LQ factorization of A
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DGELQF.335"></a><a href="dgelqf.f.html#DGELQF.1">DGELQF</a>( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
     $                INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        workspace at least M, optimally M*NB.
</span><span class="comment">*</span><span class="comment">
</span>         IF( .NOT.TPSD ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           underdetermined system of equations A * X = B
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DTRTRS.346"></a><a href="dtrtrs.f.html#DTRTRS.1">DTRTRS</a>( <span class="string">'Lower'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non-unit'</span>, M, NRHS,
     $                   A, LDA, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>            IF( INFO.GT.0 ) THEN
               RETURN
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(M+1:N,1:NRHS) = 0
</span><span class="comment">*</span><span class="comment">
</span>            DO 40 J = 1, NRHS
               DO 30 I = M + 1, N
                  B( I, J ) = ZERO
   30          CONTINUE
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DORMLQ.363"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</a>( <span class="string">'Left'</span>, <span class="string">'Transpose'</span>, N, NRHS, M, A, LDA,
     $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           workspace at least NRHS, optimally NRHS*NB
</span><span class="comment">*</span><span class="comment">
</span>            SCLLEN = N
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           overdetermined system min || A' * X - B ||
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DORMLQ.377"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</a>( <span class="string">'Left'</span>, <span class="string">'No transpose'</span>, N, NRHS, M, A, LDA,
     $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           workspace at least NRHS, optimally NRHS*NB
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DTRTRS.385"></a><a href="dtrtrs.f.html#DTRTRS.1">DTRTRS</a>( <span class="string">'Lower'</span>, <span class="string">'Transpose'</span>, <span class="string">'Non-unit'</span>, M, NRHS,
     $                   A, LDA, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>            IF( INFO.GT.0 ) THEN
               RETURN
            END IF
<span class="comment">*</span><span class="comment">
</span>            SCLLEN = M
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling
</span><span class="comment">*</span><span class="comment">
</span>      IF( IASCL.EQ.1 ) THEN
         CALL <a name="DLASCL.401"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      ELSE IF( IASCL.EQ.2 ) THEN
         CALL <a name="DLASCL.404"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      END IF
      IF( IBSCL.EQ.1 ) THEN
         CALL <a name="DLASCL.408"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      ELSE IF( IBSCL.EQ.2 ) THEN
         CALL <a name="DLASCL.411"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span>   50 CONTINUE
      WORK( 1 ) = DBLE( WSIZE )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DGELS.420"></a><a href="dgels.f.html#DGELS.1">DGELS</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
